Magnetic field dependence of the critical current in stacked Josephson junctions. 
Evidence for fluxon modes in BiiSriCaCuiOs+x mesas. 
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Modulation of the critical current across layers, h{H), of stacked Josephson junctions (SJJs) as a function of 
an applied magnetic field parallel to the junction planes is studied theoretically and experimentally for different 
junction lengths and coupling parameters. It is shown that the IJfl) patterns of long SJJs are very complicated without 
periodicity in H. This is due to interaction between junctions in the stack. This, in turn, gives rise to the existence of 
multiple quasi-equilibrium Josephson fluxon modes and submodes which are different with respect to the symmetry of 
the phase and the fluxon sequence in SJJs. The critical current of long SJJs is multiple valued and is governed by 
switching between energetically close fluxon modes/submodes. Due to this, the probability distribution of the critical 
current may become wide and may consist of multiple maxima each representing a particular mode/submode. 
Experimentally, multiple branched IJ^H) patterns and multiple maxima in the 1^ probability distribution were observed 
for Bi2Sr2CaCu208+x intrinsic SJJs, which are in a good agreement with numerical simulations and support the idea of 
having different quasi-equiUbrium fluxon modes/submodes in intrinsic SJJs. 
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INTRODUCTION. 

Properties of stacked Josephson junctions (SJJs) 
are of great interest both due to possible applications and 
from the general scientific point of view. A particular 
interest to SJJs was stimulated by the discovery of high- 

superconductors (HTSC). It is well known that the 
HTSC compounds have a layered structure with 
superconducting properties mainly confined to the Cu-O 
layers. Strong experimental evidence for that was 
provided by the observation of the intrinsic Josephson 
effect [1], which is attributed to Josephson coupling 
between atomic scale superconducting layers. However, it 
is far from being totally understood. For example, the 
superconducting 'gap' obtained from measurements of 
the intrinsic Josephson effect [1,2] is about two times 
smaller than that determined from ARPES [3] and tunnel 
measurements [4]. Furthermore, unambiguous 
observations of AC as well as DC Josephson effects 
remain to be obtained. 

In this paper we study theoretically and 
experimentally the modulation of the critical current 
across layers, IJyH), of SJJs in a magnetic field parallel to 
the junction planes. Due to the fluxoid quantization the 
critical current of the Josephson junction (JJ) should 
exhibit the well known periodic 'Fraunhofer' modulation 
of /c(/f) in parallel magnetic field, see e.g. Ref. [5], with 
periodicity 

H,=^,ILA\ (1) 
where is a flux quantum, L is the junction length, and 
A* is the effective magnetic thickness [6], 

/Ci = hi-Si-SM, (2) 
which for thin layered SJJs is equal to the interlayer 
distance, s. Here 

Ai = ti + XsiCoth{di/ Xsi) + Xsi-v\Coth[di+\l , (3) 
5, = A„cosech(if,/A„). (4) 



Here X^^ and d\ are the London penetration depth and the 
thickness of S-layers and t-^ is the tunnel barrier 
thickness. Hereafter, the subscript T of a quantity 
represents its number. Modulations of the c-axis IdH) in 
HTSC intrinsic SJJs observed so far were not well 
defined [1,7,8]. Both periodic and aperiodic modulations 
of the c-axis resistivity R{H) in YBa2Cu307.x were 
reported in Ref. [9]. However, they had a 'wrong' 
periodicity, and the authors had to assume different inter- 
and intra-unit-cell junctions with dimensions different 
from the crystal lattice space periodicity. The situation 
with low-Tc SJJs is not much better. Very complicated 
Fraunhofer patterns were observed for Nb/AlO^/Nb SJJs 
[7,10,11]. To our knowledge the only clear periodic 
modulation of the perpendicular resistivity /?(//) with 
periodicity Ho=(S>o/Ls was reported for Nb/Cu multilayers 
[12]. However, such a behavior was observed only in a 
small temperature range close to the critical temperature 
T^.. At lower temperatures the modulation in R{H) 
becomes complicated without clear periodicity and with 
pronounced hysteresis [12]. 

Theoretically, the IdH) dependence for 'small' 
SJJs, with L much less than the Josephson penetration 
depth Aj, was studied in Ref. [13]. It was shown that for 
small SJJs pure 'Fraunhofer' oscillations of IdH) similar 
to that for a single JJ should be observed with a 
periodicity given by Eq.(l). In Ref. [14] the IdH) of 
'long' SJJs, L»Aj, was studied. A monotonous power- 
law decrease of IdH) was predicted. However, this 
assumes an absolutely rigid Josephson vortex lattice. 
More realistically, the critical current in long JJs is 
determined by fluxon entrance and their redistribution in 
the junction with increasing H, see e.g. [5]. In Ref. [12] it 
was argued that complicated aperiodic IdH) patterns in 
long SJJs are caused by the asynchronous fluxon entrance 
in different junctions with the increase of H, which is a 



general property of long SJJs. In particular, fluxons in 
SJJs do not necessarily form a triangular lattice but 
rather may have more complicated arrangements as it 
can be seen from numerical simulations in Ref. [7]. 
Recently, it was shown that in SJJs there are multiple 
quasi-equilibrium fluxon configurations (modes) [6]. An 
evidence for the existence of such fluxon modes in 
Bi2Sr2CaCu208+x mesas was obtained in Ref. [15]. In the 
dynamic state the existence of multiple fluxon modes in 
SJJs would result in multiple flux-flow branches in the c- 
axis I-V curves [6] and two-dimensional collective cavity 
resonances [16]. Multiple flux-flow branches have been 
observed for Nb/Cu multilayers [12] and 
Bi2Sr2CaCu208+x mesas [17]. 

In the current paper we study fluxon modes in 
SJJs and their influence on the IJ^H) patterns. We show 
that the IJ^H) patterns in coupled long SJJs are very 
complicated and are determined by the process of 
switching between closely spaced fluxon modes, which 
occurs at the field intervals much less than Hq given by 
Eq.(l) and without clear periodicity. In addition to the 
fluxon modes defined by the number of fluxons in each 
SJJ, we found that there are certain fluxon submodes, 
which are different with respect to the symmetry of the 
phase and the fluxon sequence. The IJ^H) patterns for 
different junction lengths and coupling strengths are 
simulated numerically. The periodicity of IJ^H) is 
restored when either L becomes comparable to Aj or when 
coupling between junctions vanishes. Finally, we present 
the experimental c-axis /c(^ patterns for 
Bi2Sr2CaCu208+x mesa structures. We observed multiple 
branches on the /(,(//) patterns which are in good 
agreement with numerical simulations and support the 
idea of having multiple fluxon modes in intrinsic HTSC 
SJJs. 

THEORY 

The behavior of SJJs can be described by the 
coupled sine-Gordon equation [18] which for a double 
SJJ in the static case can be written as 

= ^inifO -^^sin( f,) - ^AIZ^IJl^ (5a) 

yciAi Ai Jci 

-2 // JclAj . , . Si ■ / 1 (A2-52) 

X]^(p ^= ^ sin((p^) -sin((p^) ^ ■ (5b) 

7ciAi Ai Ai Jci 

Here (pi 2 is the phase difference across junctions 

1 and 2, respectively, the 'prime' denotes the spatial 

derivative, J^i i and J^, are the critical and the bias current 

densities, respectively and Ay, is the Josephson 

penetration depth of JJ 1 , 
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The coupling strength of double SJJs can be described by 
the dimensionless coupling parameter 
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Fig. 1. Simulated JdH) dependencies for a long 
strongly coupled double SJJs with nonidentical 
junctions J^2=^Ji:\- In the insets, the numbers of fluxons 
in junctions 1 and 2 are shown as a function of H. 
Multiple closely spaced fluxon submodes can be seen. 

The thermodynamic equilibrium for a given H is 
achieved at the minimum of Gibbs free energy, which for 
a double SJJ can be written as [6] 
G[B) = 



i=l,2 



I,^Jci(i-cos{(p^)+^[B}iAi-2BfiBf2S2 + B}2A2] 

- — [B/1A1 + B/2A2J. (8) 

Here 2 is the magnetic induction of the fluxon. From 
Eq. (8) it is seen that G{B) is a bilinear form of Bfi 2 
which can be minimized in different ways for different 
relations between Bfi and Bf2. In other words, G{B) has a 
particular minimum for each particular fluxon mode 
(«i,«2), where «i 2 are the number of fluxons in junctions 
1 and 2, respectively. 

Fig.l shows simulated JJ^H) dependencies for a 
double stack with an overlap geometry consisting of long, 
L=10Aji, nonidentical junctions J^2=2Jqi with thin S- 
layers and strong coupling 5=0.5. The lengths are 
normalized to Aji, the current density is normalized to 
y^i, and the magnetic field is normalized to 
//p=<I>o/(?iAiAji). In the insets we show the number of 
fluxons in junctions 1 and 2 determined via the total 
phase shift, «i 2=(9'i,2(i)-?'i,2(0))/2;z:at the maximum bias 
current. Simulations were made in the following way. 
The magnetic field was swept either from zero to the 
maximum field or in the opposite direction. For a 
particular H, a solution of Eq.(5) at zero bias current was 
found using a subsequent iteration, followed by the bias 
current increase. The critical current was determined as a 
maximum bias current value for which the numerical 
procedure converged. Fig.l contains several runs with 
different initial conditions corresponding to different 
phase distributions in SJJs at H=Q and H=H,^.^^. From 
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Fig. 2. Four particular sub-branches of J^H) 
corresponding to submodes of the (3,0) mode for the 
same SJJs as in Fig.l. In inset a) spatial phase 
distributions for the submodes are shown. Number of 
fluxons in JJ 1 is shown in inset b). The existence of 
several metastable solutions with different symmetry is 
seen. 

Fig. 1 it is seen that the 7c(^^ pattern is very complicated. 
It consists of a number of distinct and closely spaced 
branches. The insets in Fig.l illustrate that those 
branches are characterized by certain fluxon 
configurations in the SJJs. Two particular JdH) branches 
corresponding to a single fluxon entrance into JJ 1 and JJ 
2, modes (1,0) and (0,1), respectively, are marked in 
Fig.l. From the insets in Fig.l it is seen that there are 
closely spaced submodes even for the same number of 
fluxons in the junctions. From our analysis of spatial 
phase distributions we have found two main reasons for 
the existence of multiple submodes: (i) The sine-Gordon 
equation for long JJs allows several metastable solutions 
having different symmetry with respect to the junction 
center, (ii) The same number of fluxons in SJJs can be 
arranged in different sequences. 

In Fig. 2 we reproduce four particular sub- 
branches of Jc{H) for the same SJJs as in Fig.l for the 
(3,0) mode. The number of fluxons in JJ 1 is shown in 
inset b), from which it can be seen that each of the four 
Jc{H) sub-branches corresponds to a particular (3,0) 
fluxon submode. In inset a) the spatial phase 
distributions for the four submodes are shown at H=70 at 
zero bias current. From Fig. 2 a) it is seen that in long 
SJJs there exist several solutions of the coupled sine- 
Gordon equation for the same H. The ambiguity in the 
solution of the sine-Gordon equation for long JJs was 
previously studied by Owen and Scalapino [19]. 
Recently, it was shown that several metastable solutions 
exist in a long JJ, which correspond to local minima of 
the free energy [20]. Solutions of this type are shown in 
Fig. 2 a). From the general properties of the sine-Gordon 
equation at zero bias current it can be shown that for a 
given solution there is a 'mirror' solution which is anti- 
symmetric with respect to the junction center. Such are 
solutions 2) and 3) in Fig. 2 a) 



Fig. 3. Three sub-branches JdH) corresponding to 
submodes of the (2,2) mode with different fluxon 
sequence for the same SJJs as in Fig. 1. In the insets, 
current distributions in JJ 1 and JJ 2 for different 
submodes are shown. 

(Pa){x)=-(p(y,{L-x). (9) 
Solutions 1) and 4) do not have mirror solutions since 

(p(r)ix)=-(P(r)iL-x), 

(p(4,ix)=-(p(4)iL-x)+27i:. 
Application of the bias current removes the symmetry 
between solutions in Eq.(9) and results in each solution 
having its own critical current, as shown in Fig. 2. 

The second reason for the existence of multiple 

submodes in long SJJs is inherent only to SJJs and is 

caused by fluxon interaction in different junctions. In 

Fig. 3, three sub-branches of JdH) for the same SJJs as in 

Fig. 1 are shown corresponding to different submodes of 

the (2,2) mode. Current distributions in JJ 1 and JJ 2 are 

shown in the insets for //=80. From the insets it is seen 

that the submodes are different with respect to the fluxon 

sequence in SJJs. We note, that the self energies of the 

submodes shown in Fig. 3 are nearly equal, which means 

that SJJs can be in any of such states with nearly equal 

probability. The number of fluxon modes and submodes 

increases rapidly with increasing total number of 

junctions and fluxons. From mathematical statistics it 

follows that the number of fluxon modes 

(indistinguishable fluxon sequences) in a stack with A' 

non-identical junctions with M fluxons in total is [6] 

{N+ M-l)! 

m = -, 7 , (10a) 

{N-1)!M! 

and the number of submodes (distinguishable fluxon 
sequences) is 

m,=N'^. (10b) 
The total number of fluxon submodes increases 
fiirther due to the existence of several solutions with 
different symmetry, as discussed above, see Fig. 2, which 
are not included in Eq.(lO). The number of submodes, 
and hence, the number of different Ji:{H) branches grows 
rapidly (exponentially) with H due to the increase of M in 
Eq.(lOb). As a result the JdH) patterns of long SJJs 
become very complicated and are determined by 
switching between closely spaced sub-branches 
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Fig.4. Jc{H) pattern for a long double SJJs with 
intermediate coupling, 5=0.31. In inset a) Jc(^ for a 
long double SJJs with small coupling is shown. In inset 
b) JciH) of strongly coupled small SJJs is shown. 

corresponding to different submodes. Such switching 
occurs at field intervals much less than the period of the 
Fraunhofer pattern in a single JJ, Eq.(l). Fig.l illustrates 
that the maximum of 7c(^ in long SJJs decrease nearly 
monotonously with H and do not have apparent 
periodicity. 

In Fig.4, the JdH) pattern is shown for a long 
double SJJ with intermediate coupling, 5=0.31. It is seen 
that multiple fluxon submodes still exist. However, the 
difference between submodes is reduced with respect to 
the strong coupling case in Fig. 1 and some weak 
periodic modulation of JdH) can be seen. With further 
decrease of coupling the difference between submodes 
disappears due to the vanishing fluxon interaction in 
neighboring SJJs. This is illustrated in inset a) of Fig. 4, 
where the JdH) pattern for a long double SJJ with thick 
S-layers and small coupling is shown. In this case the 
Jc{H) is simply defined by the minimum critical current 
of individual JJs in the stack and JdH) has a periodic 
Fraunhofer modulation. 

Periodic Fraunhofer patterns JdH) can also be 
observed for small SJJs with L comparable or less than Xj 
for arbitrary coupling strength [13]. In this case there are 
no fluxons and consequently no fluxon submodes that 
might complicate the 7c(70 pattern. This is illustrated in 
the inset b) of Fig. 4, where the JdH) pattern of strongly 
coupled small SJJs is shown. Apparent periodicity given 
by Eq.(l) is seen. A few sub-branches are still visible at 
low magnetic fields due to the finite length of the stack, 
L=2. 

EXPERIMENT 

For the study of intrinsic Josephson effect in 
HTSC, mesas with lengths 20 and 50 |a,m were fabricated 
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Fig.5. The IdH) pattern for the 20 )J,m long 
Bi2Sr2CaCu208+x mesa at r=4.2K. Circles correspond to 
increasing magnetic field and diamonds - to decreasing. 
Large symbols represent the main maxima while small 
symbols - the secondary peaks in the probability 
distributions. Multiple branches of IdH) can be seen in 
Fig.5. In the insets the probability distributions of the 
critical current are shown for a) H = 1.5 mT and h) H = 
0.45 T. 

on surfaces of Bi2Sr2CaCu20g+x single crystals [21]. The 
mesas typically contain 4-10 intrinsic SJJs. Several 
contacts on the top of the mesas allowed us to make four- 
probe measurements of I-V curves. The bias current was 
swept with a low frequency and the c-axis critical 
current, /(., was measured with a 1 )iV threshold criterion 
[22]. The magnetic field parallel to the intrinsic SJJs {ab 
plane) was generated by a superconducting coil. For each 
magnetic field value, 2048 switching events from zero 
voltage state were measured in series and were plotted as 
a probability distribution histogram, P{Ic), defining the 
probability of having a particular critical current value. 
In the insets of Fig.5 the probability distributions of the 
critical current are shown for a) H = 1.5 mT and h) H = 
0.45 T and at temperature T = 4.2 K. It is seen, that at 
low magnetic fields the probability distribution is very 
wide (0.5 mA or more). Such a wide distribution cannot 
be simply explained by thermal fluctuations or external 
noise in a single JJ. The experimental setup has been 
tested by measuring low-T;. Nb-AlO^-Nb tunnel junctions, 
which show that external noise has a negligible effect. 
This can be seen also from the comparison with a high 
field case, see inset b) in Fig.5, for which the width of the 
probability distribution is much smaller (about 2.5|a,A). 
Note, that the noise level does not decrease with 
increasing field. From inset a) in Fig.5 it is seen that the 
critical current probability distribution consists of 
multiple superimposed maxima (a certain number of 
switching events with larger /(. can be also seen in the 
inset b)). Such a behavior is consistent with the idea of 
the existence of multiple quasi-equilibrium fluxon modes 
in intrinsic SJJs with different critical currents. In this 
case each maximum corresponds to a particular fluxon 
mode. The state of SJJs is not uniquely defined by 
external conditions, but rather, it can be described only 
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statistically with a certain probability of being in any of 
the quasi-equilibrium fluxon modes, defined by the 
energy of the mode. In Fig. 5, we plotted the values of c- 
axis Ic corresponding to maxima in the probability 
distribution versus the applied magnetic field for the 20 
jO-m long mesa. The length of the mesa is much larger 
than the Josephson penetration depth which is estimated 
to be less than 1 jim. Circles correspond to increasing 
magnetic field and diamonds - to decreasing. Large 
symbols represent the main maxima while the small 
symbols - the secondary peaks in the probability 
distributions. From Fig. 5 it is seen that the IdH) pattern 
of Bi2Sr2CaCu208+x mesa is complicated and does not 
exhibit clear periodicity. However, several distinct 
branches of IdH) can be seen in Fig. 5. Similar behavior 
of IdH) has been observed for all studied mesas. The 
IdH) patterns exhibit a certain hysteresis and a prehistory 
dependence and depend on the sample geometry. All this 
is in a good qualitiative agreement with our numerical 
simulation and supports the idea of existence of multiple 
quasi-equilibrium fluxon modes in HTSC intrinsic SJJs. 

In conclusion, we have shown that the IdH) 
patterns of Bi2Sr2CaCu208+x mesas are consistent with 
the existence of Josephson effect between atomic scale 
superconducting layers so that the single crystal behaves 
as intrinsic SJJs. On the other hand, when the length of 
SJJs is larger than the Josephson penetration depth the 
behavior of IdH) is much more complicated than that of 
a single JJs and does not exhibit the periodicity in H. 
This is due to the existence of multiple quasi-equilibrium 
fluxon modes with different number of fluxons in SJJs. 
From the theoretical analysis we have found that there 
exist various fluxon submodes different with respect to 
the symmetry of the phase and the fluxon sequence. The 
critical current of long SJJs is multiple valued and is 
governed by switching between quasi-equilibrium fluxon 
modes and submodes. Due to the fact that energetically 
close modes/submodes can have different the 
probability distribution of the critical current becomes 
very wide and consists of multiple maxima, each 
representing a particular mode/submode. This is exactly 
what we have observed experimentally for long 
Bi2Sr2CaCu208+x mesas. 
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